#include "tutorialVariationSystematics.hxx"
#include "AnalysisUtils.hxx"

const bool debug = false;

/*
  A variation systematic is run before the selection, modifying the input data. The selection will then proceed 
  on that modified data, such that we can test the effect of any source variation on the final number of selected events and on 
  the distributions of any observable. As advanced, a variation systematic could have two effects:
  - Change the number of selected events because a variable we can on (or depends on it) is varied.
  - Change the value of an observable. 

  In this particular example we scale the momentum of all tracks, this will have a small effect on the number of selected events 
  (through the PID cut, which depends on momentum), because there is no direct cut on momentum, but a large effect on the distribution of 
  the muon candidate momentum

*/


//********************************************************************
tutorialVariationSystematics::tutorialVariationSystematics():SystematicVariationBase(1),SystematicSource("MomentumScale",k1D_SYMMETRIC){
//********************************************************************

  /* A variation systematic inherits from two base classes:
     - SystematicVariationBase: The 1 in the constructor is the number of systematic parameters (one in this case).
     - SystematicSource: Allows dealing with the source parameters to be propagated (reads the .dat file, etc). In this particular case the 
       .dat file will be MomentumScale.dat under psycheSystematics/vXrY/data. The name of the systematic data file is the one given to the 
       SystematicSource constructor (above). The second argument in the constructor is the type of systematics source. Thse are the available types:

       For standard weight systematics: 
       k1D_SYMMETRIC:       the systematic depends on a single observable (i.e. momentum) and its error is symmetric
       k2D_SYMMETRIC:       the systematic depends on two observables (i.e. momentum and angle) and its error is symmetric
       k1D_ASSYMMETRIC:     the systematic depends on a single observable (i.e. momentum) and its error is assymmetric
       k2D_ASSYMMETRIC:     the systematic depends on two observables (i.e. momentum and angle) and its error is assymmetric

       For efficiency-like weight systematics:
       k1D_EFF_SYMMETRIC:   the systematic depends on a single observable (i.e. momentum) and its error is symmetric
       k2D_EFF_SYMMETRIC,   the systematic depends on two observables (i.e. momentum and angle) and its error is symmetric
       k1D_EFF_ASSYMMETRIC: the systematic depends on a single observable (i.e. momentum) and its error is assymmetric
       k2D_EFF_ASSYMMETRIC: the systematic depends on two observables (i.e. momentum and angle) and its error is assymmetric

       Efficiency-like weight systematics need the values for data and MC control samples (with errors) and also the nominal for the MC analysis sample
   */

  // Read the parameters for bin 0. This particular systematic has only one bin, becaouse we consider the same 
  // momentum scale systematic for all phase space (momentum, angle, etc)
  // Since this is a 1D sysmetric systematisc source only two parameters should be read for each bin: 
  // the mean (_scale) and its error (_scaleError).
  GetParametersForBin(0, _scale, _scaleError);
}

//********************************************************************
void tutorialVariationSystematics::Apply(const ToyVariations& toy, SystBoxB& boxbase, AnaEventB& event){
//********************************************************************

  /*
    This method takes three arguments:
    - ToyVariation: Contains the variation to be applied (in number of sigmas)
    - SystBoxB: Sometimes the information needed to compute the systematic is the same for all toys (in this case all tracks with a TPC segment 
       and starting in the FGD1 FV). The operation of getting that set of tracks is time comsuming. In order not to repeat that opperation for each toy
       it is better to do it only once for each event. This information is stored in this systematics box

    - AnaEventB: Sometimes we need extra information not present in the systematics box
  */

  // In this case the event is not needed. This line avoids a warning
  (void)event;

  // Take the box for FGD from the basic box
  if(!boxbase.FGDSystBox){std::cout<<" no FGDSYST Box "<<std::endl;return;}
  SystBoxFGDB& box = *(boxbase.FGDSystBox);

  // Get the relevant tracks for this systematic
  AnaTrackB** tracks = box.TracksWithTPCInFGDFV;

  // Nothing to do If there are no relevant tracks
  if (box.nTracksWithTPCInFGDFV == 0) return;

  if (debug)
    std::cout << "tutorialVariationSystematics::ApplyVariation(): " << box.nTracksWithTPCInFGDFV << std::endl;  
  
  // loop over the relevant tracks for this systematic
  for (Int_t itrk=0;itrk<box.nTracksWithTPCInFGDFV;itrk++){
    if (debug) std::cout << itrk << std::endl;

    AnaTrackB* track = tracks[itrk];    

    // Get a reference to the momentum to be varied
    Float_t& p = track->Momentum;
    
    if (debug)
      std::cout << "p0 = " << p << std::endl;

    // Apply the momentum scale factor
    p *=(1+_scale +_scaleError*toy.Variations[0]);


    if (debug) std::cout << "p = " << p << std::endl;
  }
}

//********************************************************************
bool tutorialVariationSystematics::UndoSystematic(SystBoxB& boxbase){
//********************************************************************

  /* This method is very important. When this method returns true, the whole AnaSpill class is reseted such the all possible variations 
     are undone (we get the spill as it was after corrections). This is however very time consuming since the AnaSpill class is a very big one. 

     To speed up this action we can move back to its original value only the variables that were modified. Each systematic variation will do that operation 
     for the variables modified by it.     
  */

  // Take the FGD box from the basic box
  SystBoxFGDB* box = boxbase.FGDSystBox;

  // Check that the vector of tracks exists and that is not empty
  if(!box->TracksWithTPCInFGDFV) return false;
  if(!box->nTracksWithTPCInFGDFV) return false;

  // Get the relevant tracks for this systematic
  AnaTrackB** tracks = box->TracksWithTPCInFGDFV;

  for (Int_t itrk=0;itrk<box->nTracksWithTPCInFGDFV;itrk++){
    if(!tracks[itrk]->Original) continue;
    
    // Go back to the corrected momentum
    tracks[itrk]->Momentum = tracks[itrk]->Original->Momentum;
  }

  // Don't reset the spill to corrected
  return false;
}
